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Using both analytical and simulational methods, we study low-temperature nucleation rates in 
kinetic Ising lattice-gas models that evolve under two different Arrhenius dynamics that interpose 
between the Ising states a transition state representing a local energy barrier. The two dynamics 
are the transition-state approximation [T. Ala-Nissila, J. Kjoll, and S. C. Ying, Phys. Rev. B 46, 
846 (1992)] and the one-step dynamic [H. C. Kang and W. H. Weinberg, J. Chem. Phys. 90, 2824 
. (1989)]. Even though they both obey detailed balance and are here applied to a situation that 

does not conserve the order parameter, we find significant differences between the nucleation rates 
observed with the two dynamics, and between them and the standard Glauber dynamic [R. J. 
Glauber, J. Math. Phys. 4, 294 (1963)], which does not contain transition states. Our results show 
that great care must be exercised when devising kinetic Monte Carlo transition rates for specific 
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■ I. INTRODUCTION 

a 

o , 

' Nucleation phenomena are crucially important in many scientific disciplines, ranging from biochemistry* to earth 
sciences^ astrophysics^ and cosmology;^ and they have been studied by kinetic Monte Carlo simulations in 
electrochemistryjSi&ii materials science jifliiLiSi magnetismfiii 4 * and atmospheric science^ to mention just a few. 
However, many questions in nucleation theory are not yet resolved, and recently there has been much interest 

■ in using kinetic Ising or lattice-gas systems as theoretical and computational models of nucleation. In particular, 
(f) ' over the last decade much work has been done on their dynamical behavior at very low temperatures, both in the 

computational 16 i 17 i 18 i 19 i 20 i 21 i 22 i 23 i 24 and mathematical 25 i 26 i 27 i 28 i 29 i 30 i 31 i 32 i 33 i 34 physics communities. While this body 
of literature is divided between the Ising spin and lattice-gas formulations, we here choose to use the more symmetric 
Ising language. However, we remind the reader that the Ising spin up and down states correspond to lattice-gas 
filled and empty sites, and that the magnetic field in the Ising model is proportional to a chemical potential in the 
lattice-gas case. The explicit mappings are straightforward and can be found in Endnote l35l 

In a typical numerical nucleation experiment, the system is prepared in a metastable state with all spins antiparallel 
to the applied field, and is then allowed to evolve under a stochastic dynamic until the magnetization reaches zero on its 
way to the stable phase consisting of spins mostly aligned with the field. In the magnetic (ferroelectric) interpretation, 
this situation corresponds to magnetization (polarization) switching in an applied fieldfi 3 ^^ while in the lattice-gas 
interpretation it would correspond to submonolayer adsorption in the limit that lateral adsorbate diffusion can be 
ignoredi 36 : 37 

At very low temperatures, the dynamical behavior is influenced by the discreteness of the lattice, and it is possible 
to calculate exactly both the critical droplet (the saddle-point configuration) and the most probable path the system 
follows through configuration space during a nucleation event. The average nucleation time, defined as the inverse 
of the nucleation rate per unit system volume, has, for a two-dimensional system evolving under the Metropolis 
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d ' dynamic! 3 ^ 9 , in the low- temperature limit, been rigorously shown to be 28 

(r) = AeF. (1) 

Here, A is a nonexponcntial prefactor, T is related to the energy barrier separating the metastable state from the 
stable state, and (3 — X/k-gT where T is the temperature and fee is Boltzmann's constant (hereafter set equal to one). 
With this definition, (r) is independent of the system size. 

In kinetic Monte Carlo simulations, time is usually measured in Monte Carlo steps per site (MCSS), corresponding 
to one attempted update per site in the system. Since thermal fluctuations occur in parallel all over the system, this 
is the measure of Monte Carlo time that would be proportional to physical time. At the very low temperatures and 



relatively small system sizes considered here, the average waiting time for a nucleation event is much longer than 
the time it would take a supercritical droplet to grow to a size comparable to the size of the system. As a result, 
the change of the order parameter from its metastable to its stable value is accomplished by one single droplet - 
the first one to nucleate after the simulation is started. Since nucleation events can occur independently all over 
the system, the metastable lifetime in this single-droplet regime, operationally defined as the time until the system 
magnetization reaches zero and measured in MCSS, is inversely proportional to the system volume pitLSMML However, 
the quantity of interest here is the inverse nucleation rate (r), which is independent of the system size. To remove the 
volume dependence from the time measurements and obtain data directly applicable to (r) , we therefore here choose 
to measure simulation times in units of individual Monte Carlo steps (MCS). It is in these time units that we shall 
loosely refer to (r) as "the metastable lifetime." Lifetimes measured in the conventional MCSS units can be obtained 
through division by the number of sites in the system, N. 

Neves and Schonmann2& have shown that the critical droplet is a rectangle of overturned spins of size I X (£ — 1) 
with an extra "knob" consisting of one overturned spin on one of its long sides. The critical length I is given by 

I = [2J/\H\\ + 1 , (2) 

where [x\ denotes the integer part of x. Here, J > is the nearest-neighbor interaction constant of the Ising model, 
which will henceforth be set equal to unity, and H is the applied field in the Ising model. The critical length thus 
changes discontinuously at values of \H | such that 2/\H\ is an integer. 

It has been commonly assumed that the quantity T in Eq. Q equals the height of the energy barrier separating 
the metastable state from the equilibrium state. However, we recently showed that T actually depends on the specific 
stochastic dynamic under which the system evolves^ Only for certain dynamics, which include the Metropolis 
dynamic used in Ref. |2^, as well as the commonly used standard Glauber dynamic; 3 ^ 3 , does T actually equal the 
energy barrier i 1 ^ 42 For other dynamics that also satisfy detailed balance, and so are perfectly admissible from the 
point of view of eventually bringing the system to equilibrium, T differs from the energy barrier 4£ 

However, neither of the dynamics studied in Ref. |42j contains a microscopic energy barrier against individual spin 
flips. Since such a barrier is often needed in modeling the dynamics of physical systems, we here consider two 
specific dynamics that both possess a microscopic barrier. These are the transition dynamics approximation (TDA), 
introduced by Ala-Nissila, et al.^^*^ and the commonly used one-step dynamic (OSD)>i2*i24L2£ Both can be used 
in studies of adsorption/desorption and diffusion processes. Formal definitions of these transition rates are given in 
Sec. El 

The rest of this paper is organized as follows. The Ising lattice- gas model and the different stochastic Monte 
Carlo dynamics are introduced in Sec. [HJ and the methods used and results obtained are discussed in Sec. IIIII In 
particular, analytical calculations using a one-step Markov chain are discussed in Sec. lIII A1 computer-aided analytical 
calculations using absorbing Markov chains are discussed in Sec. IIII Bl and Monte Carlo simulations using the Monte 
Carlo with absorbing Markov chains method are discussed in Sec. IIII CI A summary and conclusions are given in 
Sec. IIVI Some preliminary analytical results on the TDA dynamic are contained in Ref. l49l 



II. MODEL AND DYNAMICS 



The square-lattice S — 1/2 Ising ferromagnet with unit interaction is defined by the Hamiltonian 

H = - ^ cr a <7p - H^a a , (3) 

where the Ising spins a a = ±1, H is the dimensionless applied magnetic field, X)(q 0) runs over all nearest-neighbor 
bonds on a square lattice, and ^ Q runs over all lattice sites. (The mapping to the lattice-gas formulation is given in 
Endnote 113) When this system develops under the standard Glauber dynamic with spin-flip rat o 39 ' 43 

Wg = 1 + exp (13 AE) 1 (4) 

where AE is the energy change that would result from the transition, T in Eq. is given bjsi&S 

r G = 8£- 2\H\(£ 2 -£+ 1) . (5) 

There is agreement in the literature that for 1 < \H\ < 2, the prefactor is given by Aq = 3/[8(£ — 1)]^ For \H\ < 1 
and such that 2/\H\ is not an integer, Ref. [23 also gives this result, while Refs. instead give Aq — 3/ [8 (I — 1) +4] . 
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The extra factor of 4 in the denominator in the latter result would be due to the effect of degenerate configurations 
with the same energy and the same number of overturned spins that occur for I > Since we in this paper perform 
numerical calculations onl y fo r \H\ > 1 (£ < 2), we have not explicitly tested this discrepancy. Results for \H\ > 2 were 
obtained in Refs. l20EHl22ll42l Note that Tq is a continuous function of \H\, even though £ and Aq are discontinuous 
wherever 2/|i?| is an integer. 

The interpretation of T as given in Eq. J^J) is simply the energy difference between the metastable state and the 
critical droplet. In Ref. |42jwe showed that there exist stochastic dynamics that satisfy detailed balance, as does the 
Glauber dynamic, but for which T is different from the value given by Eq. (J3J, even though the critical droplet remains 
unchanged by the change of dynamic. The specific dynamic for which this was demonstrated was the "soft Glauber 
dynamic,"^ for which the transition rates are given by 

WsG = 1 + exp (0AEj) 1 + exp ((3AE H ) ' ^ 

where AEj and AEh are the changes in the interaction energy and in the field energy, respectively. Consequently, 
the total energy change is AE = AEj + AEh- It is the factorization property into one part that depends only on 
AEj and another that depends only on AEh that defines this dynamic as "soft." 50,51,52 « In contrast, the standard 
Glauber dynamic, Eq. (@J, which does not have this factorization property, is described as "hard." The need to use 
soft dynamics in cases where the 'field' represents a chemical-potential difference has been recognized in some Monte 
Carlo simulations of crystal growthiSiiiS 

In many applications it is physically reasonable to introduce a microscopic transition barrier^&Siii Such a barrier 
is not part of the Ising or lattice-gas Hamiltonian, and it is therefore absent in the transition rates discussed above. 
Rather, the barrier represents a transition state which is inserted between the states allowed in the Hamiltonian, such 
as a saddle point in a corrugation potential for particle diffusion, or a high energy associated with a transitional spin 
state that is not along one of the two directions allowed by the Ising Hamiltonian. Dynamics that contain such a 
barrier are often called Arrhenius dynamics, and they can be viewed as a way to use the lattice-gas or Ising model 
to simulate dynamics in an underlying, continuous potential. 7,45 Here we use the following approximation to express 
the transition-state energy i^ 4 ^ 5 ^^ 

E^^ + A, (7) 

where Ei and E{ are the initial and final energies, and A is the bare, microscopic energy barrier. In electro- 
chemical applications, such as electron- or ion-transfer reactions, this corresponds to the symmetric Butler- Volmer 
approximationi^SiS The construction corresponding to Eq. J7J) is illustrated in Fig. ^ 

The two specific dynamics considered in this paper are the TDA and OSD rates with the transition-state energy 
E T given by Eq. 0. They are defined ao 44 i 45 i 46 

1 1 

WTDA ~ I + exp[(3(E T - E,)} I + exp[/3(£ f - E T )} ' (8) 

andi 7 -^ 

Wosd = exp[-0(E T - Ei)] = exp(-/?A) exp (-0 ® ~ , (9) 

respectively. From the definitions of hard and soft dynamics given above, we see that the TDA dynamic cannot be 
factored into parts that depend only on AEj and AEh and thus is hard. On the other hand, the OSD dynamic can be 
factored in this way and thus is soft. While these dynamics are often used for diffusion problems (i.e., conserved order 
parameter) we will here use them in the simpler context of magnetization switching or adsorption/desorption 

(i.e., nonconserved order parameter). 



III. RESULTS 

Like we did in Ref. l42l in this paper we obtain our results in three different ways. First, we calculate analytically 
by hand the first-passage time from the metastable state to an absorbing state just beyond the saddle point in an 
approximation that the path in configuration space corresponds to a simple one-step Markov process^ Second, we 
perform computer-aided analytic calculations using the technique of absorbing Markov chains (AMC)^^^ allowing 
for multiple branching paths and "blind alleys" on the way to the absorbing state(s). Third, we perform Monte Carlo 
simulations using the Monte Carlo with AMC (MCAMC) techniqueii 7 ^^ 1 ^ The first method provides the clearest 
physical insight, and for noninteger values of 2/\H | the results are fully confirmed by the other two. 
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A. One-step Markov chain 



As we discussed in Ref. [42, the one-step Markov chain for 1 < H < 2 {I = 2) corresponds to the configurations 
labeled i — 0, ...,4 in Fig. The label i corresponds to the number of overturned spins, such that the starting 
configuration with all spins up has i = 0, and the saddle point has i = i* = 3. In general the absorbing state is 
labeled I > i* > 1. The mean time spent in state i (the sojourn time) is denoted hi, the rate at which the cluster of 
overturned spins grows from i to i + 1 is denoted gi, and the rate with which it shrinks from i to i — 1 is denoted Si. 
These quantities satisfy the equation 58 ' 63 - 64-65 

h i -i = (s i h i + N)/g i -i (10) 

with boundary conditions si = sq = 0. The number N of sites in the system represents the total probability current 
through the Markov chaini 5 ^ This relation enables us to obtain hi recursively as hj^i = N/gi^i and 

9i £r[ gi+k fJi \9i+j-iJ 

for < i < I — 2. Assuming that I > i* and that the growth time for the supercritical droplet can be neglected 
in comparison with the nucleation time, the metastable lifetime (r) is given by the mean first-passage time to /, 
(t) = (tj) = J2i=o hi- Grouping the terms in the sum according to the "unpaired" factors l/gi+fc in Eq. Qllfl. we get 

fr>=-+E-(i+£ff— I ■ 

With a suitable interpretation of the factor N, this result is general for any one-step Markov chain with absorption 
at /, regardless of the values of gi and Sj. However, for the Ising or lattice-gas model the transition probabilities are 
related by detailed balance, so that 



S il n i = e fi{E -L i 



(13) 



where Ei is the energy of state i. The degeneracy factors nf and n\_-i are the numbers of lattice sites at which a 
single spin flip can shrink the cluster from i to i — 1 and analogously for growth from i — 1 to i, respectively. As a 
result, Eq. I|12fl becomes 

fa>=-+£- fi+E ^ Ei - Ei - k) n ( i4 ) 

5o tt 3i \ fx to n ^-i / 

The number of sites that can give rise to a transition from zero to one overturned spin equals the system size A, so 
rip = N. In contrast, the other nf and n| are just functions of the cluster's size (i) and its shape, and thus independent 
of N. Consequently, each term in Eq. (|14fl that has a factor of go or in the denominator is independent of N. 
This turns out to be the first term, N/go, as well as exactly one term for each value of I in the outer sum, namely, 
the one corresponding to k = ! in the inner sum. All the other terms are of O(N). We note that the A- independent 
terms in Eq. Ijl4|l constitute a series analogous to the expression for the inverse nucleation rate, given in Eq. (15) of 
Ref. |2^ In the limit of /3 — > oo, Eq. (|14f> is dominated by the term or terms with the largest exponential factor. The 
selection of these terms, which determine both T and A, is described below. However, we first need to determine the 
low-temperature spin-flip rates in the two different dynamics considered here. 

For the square-lattice Ising system, all possible spin configurations can be classified into 10 classes, determined by 
the value a of a spin (+ for the metastable direction and — for the stable direction) and the number A + of its nearest 
neighbors that are in the metastable direction^ In the low-temperature limit, the rate p m for flipping a spin in class 
m (m = 1,2,.. .,10) [given by Eqs. (JSJ) and @ for the TDA and OSD dynamics, respectively], can be simplified as 
shown in Tabled 

Figure [3 corresponds to a one-step Markov chain with I = 4. For 1 < \H\ < 2 the saddle-point configuration 
(i = i*) is the L-shaped three-site cluster with I = 2 (i* = 3). Among the dominant terms in Eq. (|14fl is always the 
term with k = I = i* — 1. For all H < 2, the growth from i* — 1 to i* always involves adding a "knob" to the long side 
of an £ x (£— 1) rectangle, such that <7j»_i = nf*_ 1 P2 = 2£p 2 - For 2 < \H\ < 4, the saddle point is a single overturned 
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spin (£ — 1), so i* — 1 and gi--\ = go = n^px = Np%. To avoid TV-dependent growth terms in Eq. 114(1. for > 4 
the lifetime must be defined as the first-passage time to one overturned spin, so that (r) = (n) = N/g = 1/pi- 

To obtain T and A in the two dynamics for \H\ > 1, we explicitly write out Eq. (|14[1 for 1 — 4 with the Si, gi-i, 
and Ei obtained from Fig. |2 

(r 4 > = l + ^fiV + e ^ 8 - 2 l ff l)U— fTV + ^e^^D + ie^ 12 - 4 ^^ 
Pi 4p 2 \ / 4p 2 V 2 2 / 

+ — (n+ — e ^- 2 \ H \) + E e W-MH\) + I e /5(16-6|H|)\ 
P3 V 2 4 4 ' / 

= A + B + C + V. (15) 

Here, „4 is the first-passage time from i = to 1, B is the first-passage time from i = 1 to 2, etc. Using the values of 
p m from Table [I] we can identify the dominant terms in (T4) and explicitly obtain A and T for both dynamics for all 
\H\ > 1. To avoid complicating overlaps of different field regimes, in the remainder of this paper we restrict ourselves 
to < A < 1. 

TDA dynamic: For all \H\ < 4, (T4} is dominated by one or more of the Af-independent terms. For 1 < \H\ < 2 — A 
the dominant part is V + C, for 2 - A < \H\ < 2 it is C, for \H\ = 2 it is C + B, for 2 < \H \ < 2 + A it is B, for 
2 + A < \H\ < 4 — A it is B + A, and for 4 - A < \H\ < 4 it is A. As mentioned above, for \H\ > 4 the system is 
unstable, so the terms B...D, whose dominant parts are then of O(N), correspond to growth of a supercritical droplet, 
rather than nucleation. Consequently, the metastable lifetime must be redefined as the first-passage time to a single 
overturned spin, (r) = (rj) = A, which is still independent of N. The corresponding results for Ttda and ^tda for 
all values of \H \ > 1 are listed in Table ITT1 and shown in Figs. |3| and 0] 

OSD dynamic: For all \H\ < 4 + A, (t^) is dominated by one or more of the ^-independent terms. For 1 < \H\ < 2 
the dominant part is C, for \H \ = 2 it is C + B, for 2 < \H\ < 3 it is B, for \H\ = 3 it is B + A, and for 3 < \H\ < 4 + A 
it is A. In the same way as for the TDA dynamic for \H\ > 4, for \H\ > 4 + A, the 0(N) growth terms in B...V 
would dominate the nucleation term A and must be removed by defining (r) = (tj) = A. The corresponding results 
for Fqsd and Aqsd for all values of \H\ > 1 are listed in Table ITU and shown in Fig.JSJ 



B. Computer-aided analytic calculation 



To check the results of the one-step Markov-chain approximation, analytic calculations of (r) were also performed 
by the AMC method^Sii using t = 13 transient states and all the a absorbing states that can be reached from the 
transient subspace by a single spin flip. The transition matrix defining the AMC has the form 

M (a+t)x(a+t) = (^;; Q £««J , (i6) 

where the sizes of the blocks have been explicitly indicated. Here, I is an identity matrix, is a matrix with all zero 
elements, T is called the transient matrix, and R is called the recurrent matrix. The latter two matrices contain the 
transition rates within the transient subspace and from the transient to the absorbing subspace, respectively. We use 
the convention (common in mathematics and statistics) that a probability density vector multiplies M from the left. 
Therefore, in order for M to be a Markov matrix, each of its row sums must equal one. Given that the system is 
initially described by a probability density represented by the row vector (here taken to correspond to probability 
one for the metastable state with no overturned spins, transient state in Fig.EI, the mean first-passage time to the 
absorbing subspace can be shown to be^SS 

(r> = ^(I-T)- 1 e, (17) 

where eis a column vector with all elements equal to one, and I is the t xt unit matrix. The transient and absorbing 
states are shown in Fig. [5] The transition rates that form the elements of T were calculated using the pi from Table UJ 
and degeneracies nf and n\ deduced from the geometric configurations shown in Fig. |SJ Equation l|17f) was then 
solved analytically for (3 — > 00 in a computer-aided calculation using Mathematical The metastable lifetime was 
calculated in this way, both at the special values of \H\ where the one-step analytical calculations indicated that A 
should have discontinuities (\H\ — 1.5, 2.0, 2.5, 3.5, and 4.5 for the TDA dynamic, and \H\ — 2.0 and 3.0 for the 
OSD dynamic), and at several points in the field intervals where A was expected to be constant. The results of these 
analytic calculations are included in Tables ITTI and ITTTI and also shown in Figs. 0] and [SJ The only field value where 
differences were found between the one-step and full analytic calculations was \H | = 2 for both dynamics. 
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C. Monte Carlo simulations 

Both sets of analytic results were further checked by continuous-time Monte Carlo simulations with time measured 
in MCS as elsewhere in this paper. For \H\ < 4.5 we used the Monte Carlo with absorbing Markov chains (MCAMC) 
method^*£2i£ii££ which significantly speeds up the simulations while faithfully preserving the underlying dynamic. 
The MCAMC method was applied with two transient states (states TO and Tl in Fig. HJ). After exit from this 
transient subspace, the simulation continued with the standard n-fold way algorithm£ii£& until the magnetization 
reached zero. At the lowest temperatures the MPFUN arbitrary-precision packag o^?i9? was used to obtain sufficient 
numerical precision in the inversion of (I — T) in Eq. I|17|). which in this case is a 2 x 2 matrix. For \H\ > 4.5, 
where the lifetime is the average first-passage time to a single overturned spin, we used the n-fold way algorithm, 
needing only to consider the transitions out of class with rate Npi . The system was an L — 24 square with periodic 
boundary conditions, and 2000 escapes were used (100 000 for \H\ > 4.5). Simulations were performed, both at the 
same special values of \H that were studied by the computer-aided analytic AMC method, and at several points in 
the field intervals where A was expected to be constant. The system was started in the metastable state with all 
spins up and subjected to a negative field of magnitude \H\, and the metastable lifetime was defined as the average 
time until the system magnetization reached zero (except for \H\ > 4.5, where it was defined as the average time 
until a single spin was overturned). The lifetime was then measured at each \H\ for a series of temperatures between 
T = 0.08 and a field-dependent minimum value (between approximately 0.04 for the smallest \H\ and about 0.003 
for the largest \H\). Some of the results, shown vs T as Tln(r) Y + Tin A, are plotted in Fig. for the TDA 
dynamic. To find T and A, linear fits to such plots were performed, deleting higher- T points (which may deviate 
from the low-temperature linear T-dependence) until the Q-factor of the fit became acceptable^ The results of these 
simulations are included in Tables ITT1 and 11111 and also shown in Figs. 0] and [SJ The agreement with the analytical 
results is everywhere excellent. At \H\ = 2, where there are differences between the one-step and full analytical results 
for A, the simulated results agree with the full analytical ones to within one standard deviation for both dynamics. 

IV. SUMMARY AND CONCLUSIONS 

In this work we have continued our studies of the influence of different kinetic Monte Carlo transition rates 
on the simulated structure and mobility of interfaces£2i£ii£iiZ£ and on the low-temperature simulated nucleation 
propertiesi2i2i2ii^ of Ising and lattice-gas models. These studies have documented that different transition rates can 
lead to very significant differences in both structure and dynamics, even between dynamics that all satisfy detailed 
balance and also have the same conservation properties. As an extension of these studies, we here considered two 
popular dynamics that include a local energy barrier representing a transition state inserted between individual Ising 
or lattice-gas states. Such Arrhenius dynamics, as they are often called, are appropriate in kinetic Monte Carlo simu- 
lations of discrete Ising or lattice-gas models in which the discrete states serve as approximations for high-probability 
configurations in an underlying continuous potential. Examples are the use of a lattice-gas model to study diffusion 
in a continuous corrugation potential^^^& or the use of an Ising model as an approximation for a continuum spin 
model with a strong uniaxial anisotropy^ 

The two Arrhenius dynamics that we studied are the commonly used one-step dynamic (OSD) 1 ^ 4 ^^ [Eq. J5J] and 
the two-step transition dynamics approximation (TDA) 44-45 46 [Eq. (JSJ)]. Using the same analytical and numerical 
methods that we used in Ref . 1421 for the standard Glauber dynamic [Eq. J3J] and the soft Glauber dynamic [Eq. ©], 
we calculated the exponential argument T and prefactor A in Eq. (JTJ for the metastable lifetime (r) (the inverse of 
the nucleation rate per unit volume) in nucleation simulations at very low temperatures. (Although the dynamics 
studied are different, the more detailed descriptions of the methods given here can also serve as supplementary 
material for the necessarily brief descriptions in Ref. ^3-) While the differences between the OSD and TDA dynamics, 
and between them and the standard Glauber dynamic, are not as dramatic as those between the standard and soft 
Glauber dynamics in Ref.l42l they are definitely significant. Once again they confirm that, contrary to what has been 
a common part of the folklore of nucleation theory, T does not always equal the height of the energy barrier separating 
the metastable state from the stable state. The differences that we find between the results for the OSD and TDA 
rates also indicate that these dynamics are not physically equivalent. It is thus interesting to note that differences 
have also been observed between the TDA and a version of the OSD in which Et is constant (called the initial-value 
approximation, IVA) in simulations of diffusion in a lattice-gas model of oxygen on the (100) surface of tungsten at 
elevated temperatures i^SiffiiZLZSiS 

One way of obtaining the correct dynamics for a given system could be to start from the quantum-mechanical 
formulation, and from that derive the dynamic for the classical model within the appropriate approximations relevant 
for the particular system. For example, one can utilize the time-dependent quantum density matrix for a Hamiltonian 
coupled to a heat bath, integrate out the heat bath, and with additional approximations usually arrive at a dynamic 
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for a classical or semiclassical system. This was accomplished for the quantum spin-1/2 system with a fermion heat- 
bath, giving the Ising model with the standard, hard Glauber dynamic^ It has recently also been accomplished for 
a quantum spin-1/2 lattice system coupled to a bosonic heat bath, giving a dynamic that is different from the hard 
Glauber dynamici 19 i 23 i 24 

In this paper we have only studied the effect that changing the dynamic has on the two-dimensional square-lattice 
Ising model. It is natural to ask how widely such influences on nuclcation occur in other models and in other 
dimensions. At low temperatures, studies have also been made of the three-dimensional Ising model, both through 
simulations^ an d through rigorous mathematical analysisiS^ii^ 2 ^ (Note that the res ults gi ven for three-dimensional 
systems in Refs. I25l82l83l are not correct for all \H\: correct results are given in Refs. Il8l8l! 1 We expect that both T 
and the prefactor should also change for different dynamics for the cubic-lattice Ising model. In fact, we anticipate 
that results similar to those presented here may hold for most lattice models with discrete state spaces. One important 
question is whether in off-lattice simulations, for example ones that use energy barriers from zero-temperature first- 
principles calculations^ a change in the dynamic could also change both T and the prefactor in low-temperature 
nucleation. If this is the case, then, since a small change in T can lead to orders-of-magnitude changes in lifetimes, 
particular attention must be paid to deriving the dynamic to be used in off-lattice simulations. 

The results of this study reinforce the message that great caution must be exercised when constructing transition 
rates for kinetic Monte Carlo simulations of specific physical systems 42i5L5iiIiiZ2i On the other hand, the significant 
differences between individual dynamics that have been observed indicate that experimental results for such quantities 
as nucleation rates and interface mobilities should provide valuable input for devising physically sound kinetic Monte 
Carlo algorithms in the future. 



Acknowledgments 



G.M.B. gratefully acknowledges the hospitality of the Mississippi State University Department of Physics and 
Astronomy and ERC Center for Computational Sciences, and of the Florida State University School of Computational 
Science and Information Technology and Center for Materials Research and Technology. This work was supported by 
NSF grants No. DMR-0120310 and DMR-0240078. 



* Electronic address 
t Electronic address 

* Electronic address 
§ Electronic address 



17 
18 

1!) 
20 
21 



Theory, Experiment, and 



lDuendia@ush.-gj] 
rikvold@csit.fsu.edu 
park@dave.nrl.navy.mil 
man40@ra.msstate.edu 

J. Onuchic, Z. Luthey-Schulten, and P. G. Wolynes, Annu. Rev. Phys. Chem. 48, 545 (1997). 
A. C. Lasaga, Kinetic theory in the earth sciences (Princeton Univ. Press, Princeton, 1998). 
A. B. C. Patzer, A. Gauger, and E. Sedlmayr, Astron. Astrophys. 337, 847 (1998). 
H. A. Kastrup, Phys. Lett. B 419, 40 (1998). 
H. A. Kastrup, Ann. Phys. Berlin 9, 503 (2000). 

G. Brown, P. A. Rikvold, S. J. Mitchell, and M. A. Novotny, in Interfacial Electrochemistry: 
Application, edited by A. Wieckowski (Marcel Dekker, New York, 1999), p. 47. 

S. J. Mitchell, S. Wang, and P. A. Rikvold, Faraday Discuss. 121, 53 (2002). 

F. Berthier, B. Legrand, J.Creuze, and R. Tetot, J. Electroanal. Chem. 561, 37 (2004). 

F. Berthier, B. Legrand, J.Creuze, and R. Tetot, J. Electroanal. Chem. 562, 127 (2004). 

N. Combe, P. Jensen, and A. Pimpinelli, Phys. Rev. Lett. 85, 110 (2000). 

S. Auer and D. Frenkel, Nature 409, 1020 (2001). 

K. A. Fichthorn, M. L. Merrick, and M. Schemer, Appl. Phys. A 75, 17 (2002). 
M. A. Novotny, G. Brown, and P. A. Rikvold, J. Appl. Phys. 91, 6908 (2002). 

H. L. Richards, S. W. Sides, M. A. Novotny, and P. A. Rikvold, J. Magn. Magn. Mater. 150, 37 (1995). 
R. Mahnke, R. Kaupuzs, and V. Frishfelds, Atmos. Res. 65, 261 (2003). 

M. A. Novotny, in Computer Simulation Studies in Condensed Matter Physics IX, edited by D. P. Landau, K. K. Mon, and 

H.-B. Schiittler ( Springer- Verlag, Berlin, 1997), p. 182. 

M. A. Novotny, Computer Phys. Commun. 147, 132 (2002). 

M. A. Novotny, in Computer Simulation Studies in Condensed Matter Physics XV, edited by D. P. Landau, S. P. Lewis, and 
H.-B. Schiittler (Springer- Verlag, Berlin, 2003), p. 7. 

K. Park, M. A. Novotny, and P. A. Rikvold, Phys. Rev. E 66, 056101 (2002). 
V. A. Shneidman, J. Stat. Phys. 112, 293 (2003). 

V. A. Shneidman and G. M. Nita, Phys. Rev. Lett. 89, 025701 (2002). 
V. A. Shneidman and G. M. Nita, Phys. Rev. E 68, 021605 (2003). 



7 



23 K. Park and M. A. Novotny, Computer Phys. Commun. 147, 737 (2002). 

24 K. Park and M. A. Novotny, in Computer Simulation Studies in Condensed Matter Physics XIV, edited by D. P. Landau, 
S. P. Lewis, and H.-B. Schiittler (Springer- Verlag, Berlin, 2002), p. 134. 

25 A. Bovier and F. Manzo, J. Stat. Phys. 107, 757 (2002). 

26 P. Dehghanpour and R. H. Schonmann, Commun. Math. Phys. 188, 89 (1997). 

27 P. Dehghanpour and R. H. Schonmann, Prob. Theory and Related Fields 107, 123 (1997). 

28 E. Jordao Neves and R. H. Schonmann, Commun. Math. Phys. 137, 209 (1991). 

29 R. Kotecky and E. Olivieri, J. Stat. Phys. 75, 409 (1994). 

30 E. Olivieri and E. Scoppola, J. Stat. Phys. 79, 613 (1995). 
R. H. Schonmann, Commun. Math. Phys. 147, 231 (1992). 
R. H. Schonmann, Commun. Math. Phys. 161, 1 (1994). 
E. Scoppola, J. Stat. Phys. 73, 83 (1993). 

E. Scoppola, Physica A 194, 271 (1993), and references cited therein. 

Starting from the Ising Hamiltonian of Eq. JHJ, the standard, explicit mapping between the Ising and lattice-gas formulations 
is as follows. We identify the Ising variable a a = +1 ("spin up") at site a with the lattice-gas variable c a = 1 (occupied or 
"solid") and a a — — 1 with c a = (empty or "fluid"), so that a a — 2c a — 1. The Ising and lattice-gas interaction constants 
[J and <f) respectively; J is set equal to unity in Eq. J3J] are related as (j> = 4J, and the applied Ising field H is related to 
the lattice-gas chemical potential \i as H = (/i — (M>)/2, where fio = ~ 8 J = ~2(j> is the coexistence value of \i. 

36 R. A. Ramos, P. A. Rikvold, and M. A. Novotny, Phys. Rev. B 59, 9053 (1999). 

37 M. A. Novotny, P. A. Rikvold, M. Kolesik, D. M. Townsley, and R. A. Ramos, J. Noncryst. Solids 274, 356 (2000). 

38 N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953). 

39 D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University Press, 
Cambridge, 2000). 

40 P. A. Rikvold and B. M. Gorman, in Annual Reviews of Computational Physics I, edited by D. Stauffer (World Scientific, 
Singapore, 1994), p. 149. 

41 P. A. Rikvold, H. Tomita, S. Miyashita, and S. W. Sides, Phys. Rev. E 49, 5080 (1994). 

42 K. Park, P. A. Rikvold, G. M. Buendfa, and M. A. Novotny, Phys. Rev. Lett. 92, 015701 (2004). 

43 R. J. Glauber, J. Math. Phys. 4, 294 (1963). 

44 T. Ala-Nissila, J. Kjoll, and S. C. Ying, Phys. Rev. B 46, 846 (1992). 

45 T. Ala-Nissila and S. C. Ying, Prog. Surf. Sci. 39, 227 (1992). 
" T. Ala-Nissila, R. Ferrando, and S. C. Ying, Adv. Phys. 51, 949 (2002). 

H. C. Kang and W. H. Weinberg, J. Chem. Phys. 90, 2824 (1989). 

48 K. A. Fichthorn and W. H. Weinberg, J. Chem. Phys. 95, 1090 (1991). 

49 G. M. Buendfa, P. A. Rikvold, K. Park, and M. A. Novotny, Rev. Mex. Ffs. (in press). 
P. A. Rikvold and M. Kolesik, J. Phys. A 35, L117 (2002). 
P. A. Rikvold and M. Kolesik, J. Stat. Phys. 100, 377 (2000). 

J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, 
1999). 

H. Guo, B. Grossmann, and M. Grant, Phys. Rev. Lett. 64, 1262 (1990). 
54 M. Kotrla and A. C. Levi, J. Stat. Phys. 64, 579 (1991). 

F. Hontinfinde, J. Krug, and M. Touzani, Physica A 237, 363 (1997). 
W. Schmickler, Interfacial Electrochemistry (Oxford Univ. Press, New York, 1996). 
J. O. Bockris and A. K. N. Reddy, Modern Electrochemistry, Vol. 2 (Plenum, New York, 1970). 

N. G. Van Kampen, Stochastic Processes in Physics and Chemistry, 2nd. Ed. (North-Holland, Amsterdam, 1992), Chap. VI. 
M. Iosifescu, Finite Markov Processes and their Application (Wiley, New York, 1980), p. 99. 

M. A. Novotny, in Annual Reviews of Computational Physics IX, edited by D. Stauffer (World Scientific, Singapore, 2001), 
p. 153. 

M. A. Novotny, Phys. Rev. Lett. 74, 1 (1995); erratum 75, 1424 (1995). 

M. A. Novotny and S. M. Wheeler, in Computer simulations of surfaces and interfaces, edited by B. Diinweg, et al. (Kluwer 
Academic Publishers, Amsterdam, 2003), p. 225. 

M. Kolesik, M. A. Novotny, P. A. Rikvold, and D. M. Townsley, in Computer Simulation Studies in Condensed Matter- 
Physics X, edited by D. P. Landau, K. K. Mon, and H.-B. Schiittler (Springer- Verlag, Berlin, 1998), p. 246. 
M. Kolesik, M. A. Novotny, and P. A. Rikvold, Phys. Rev. Lett. 80, 3384 (1998). 
M. A. Novotny, M. Kolesik, and P. A. Rikvold, Computer Phys. Commun. 121-122, 330 (1999). 
A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, J. Comput. Phys. 17, 10 (1975). 
S. Wolfram, The Mathematica Book, 3rd Ed. (Cambridge Univ. Press, Cambridge, 1996). 
D. H. Bailey, NASA technical report RNR-90-022. 
D. H. Bailey, ACM Trans. Math. Software 21, 379 (1995). 

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran, Second Ed. (Cambridge 
Univ. Press, Cambridge, 1992). 

P. A. Rikvold and M. Kolesik, Phys. Rev. E 66, 066116 (2002). 
P. A. Rikvold and M. Kolesik, Phys. Rev. E 67, 066113 (2003). 
J. D. Munos, M. A. Novotny, and S. J. Mitchell, Phys. Rev. E 67, 026101 (2003). 

I. Vattulainen, J. Merikoski, T. Ala-Nissila, and S. C. Ying, Phys. Rev. Lett. 79, 257 (1997). 



47 



8 



C. Uebing and V. P. Zhdanov, Phys. Rev. Lett. 80, 5455 (1998). 

I. Vattulainen, J. Merikoski, T. Ala-Nissila, and S. C. Ying, Phys. Rev. Lett. 80, 5456 (1998). 
C. Uebing and V. P. Zhdanov, J. Chem. Phys. 109, 3197 (1998). 

I. Vattulainen, S. C. Ying, T. Ala-Nissila, and J. Merikoski, J. Chem. Phys. Ill, 11232 (1999). 

C. Uebing and V. P. Zhdanov, J. Chem. Phys. Ill, 11234 (1999). 
P. A. Martin, J. Stat. Phys. 16, 149 (1977). 

G. Ben Arous and R. Cerf, Electronic J. Probability 1, paper No. 10, 1 (1996). 

D. Chen, J. Feng, M. Qizn, Science in China (Series A) 40, 832 (1997). 
D. Chen, J. Feng, M. Qizn, Science in China (Series A) 40, 1129 (1997). 



9 



TABLE I: Rates of flipping a spin in class m, p m , for \H\ > and A > in the limit of (5 — > oo. Here a = + (— ) corresponds 
to a spin in the metastable (stable) direction, and JV+ is the number of its nearest-neighbor spins in the metastable direction. 
The analytical forms of the OSD results do not depend on \H\. The TDA results marked * all take the value \ for A = 0. The 
first two lines of TDA results for classes 4 and 9 are only applicable for A > 2, and the first two lines of TDA results for classes 
5 and 10 are only applicable for A > 4. 
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TABLE II: Analytical and simulated j4tda and Ttda for the TDA dynamic with \H\ > 1. Theoretical results are valid for 
< A < 1, while all numerical values are for A = 0.5. For Atoa, the analytical AMC calculation with a one-step Markov 
chain and the full AMC calculation may give different results for fields H such that 2/\H\ is a positive integer, as well as for 
\H\ ± A. They are therefore given as ^4TDA(l-step) and Atda^uII), respectively. For Ttda, the two analytical methods always 
agree, and the results are given together as FxDA(anaL). The simulated results, ArDA(sim.) and rTDA(sim.), are obtained 
from weighted two-parameter fits to plots of Tln(r) ~ r + Tin A, such as those shown in Fig. Standard errors in the last 
digit are given in parentheses. 
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TABLE III: Analytical and simulated Aosd and Tosd for the OSD dynamic with \H\ > 1. Theoretical results are valid for 
< A < 1, while all numerical values are for A = 0.5. The estimation procedures and the meanings of all symbols are the 
same as in Table IH1 
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initial transition final 
state state state 

FIG. 1: Schematic picture of the transition barrier in the symmetric Butler- Volmer approximation, used to calculate the TDA 
and OSD transition rates. 

£-0 E 2 =12-4\H\ E 4 =16-8\H\ 

g Q =Np { g { = 4p 2 g 2 =4p 2 g 3 = p 3 




s i = Pb s 2= 2 Pi s 3 =2p 7 
E [ =8-2\H\ E 3 =16-6\H\ 



FIG. 2: The states in the one-step Markov chain of clusters of i = 0, ...,4 overturned spins, used to calculate the metastable 
lifetime (r) analytically by hand. The right-pointing arrows give the growth rates gi-i, and the left-pointing arrows give the 
shrinkage rates Si for i = 1, 2, and 3. The numbers preceding the spin-flip rates p m (see Table are the degeneracy factors 
n l-i an d n t for 9i-i an d Si, respectively. The energies Ei (relative to the metastable state, i = 0) are given at the top of the 
figure for even i and at the bottom for odd i. See discussion in the text. From Ref. l42l 
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FIG. 3: Analytical (lines) and simulated (data points) results for Tln(r) « T + Tin A for the TDA dynamic, shown vs T for 
three different values of \H\. Error bars, of the order of T/V2000, are too small to be seen on the scale of this figure. Weighted 
two-parameter linear fits to the simulation data yield the simulated estimates for A and T, given in Table ITT1 and shown in 
Fig.H 




FIG. 4: Analytical and simulated results for V (a) and A (b) vs \H\ for the TDA dynamic, using A = 0.5. Results for F for 
the standard Glauber dynamic (from Ref. 4^) are shown by dashed lines in (a). See discussion in the text. 
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FIG. 6: The transient (T) and absorbing (A) states used in the full AMC calculations with 13 T and 13 A states. Squares denote 
sites with an overturned spin. For states that do not consist of a contiguous cluster of nearest-neighbor sites, connecting dashed 
lines indicate that the overturned spins are separated by a single site [T state 4 (T4) and A state 1 (Al)]. In the noncontiguous 
state T5, the two overturned spins can be at any sites that are further apart than third-nearest neighbors. Among the A 
states, noncontiguous sites that are not connected by dashed lines can be in any relative positions, except nearest-neighbor. 
Several of the A states can also have other shapes than the ones shown here. For example, the four-site cluster in A7 can also 
be L-shaped like T10, and the low-symmetry five-site clusters A6 and A9 can have several different shapes. In addition, all 
states also include configurations obtained from those shown by rotation or reflection. All the A states, including variations 
and states equivalent under symmetry, were taken into account in the calculation. 
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